library(sf)
library(raster)
library(osmdata)
library(rmapshaper)

# Perimètre d'analyse autour du centre
buf <- 30000

# Centres CNRS
cent <- read.csv(file = "data-raw/cnrs.csv")

# Géométries LAU et population
com <- st_read(dsn = "data-raw/LAU_RG_01M_2019_4326.shp", quiet = TRUE)
country <- aggregate(x = com[,"CNTR_CODE"], by = list(com$CNTR_CODE),
                 FUN = head, 1)
com <- st_transform(com, 2154)
country <- st_transform(country, 2154)

# Sélectionner la France et généraliser ses contours
fr <- country[country$CNTR_CODE == "FR",]
fr <- st_cast(fr, "POLYGON")
fr$id <- row.names(fr)
fr$area <- st_area(fr)
fr <- fr[order(fr$area, decreasing = TRUE),]
fr <- fr[1, "CNTR_CODE"]
fr <- ms_simplify(fr, keep = 0.05)

write_sf(obj = fr, dsn = paste0("data/fr.gpkg"),
           delete_layer = TRUE, append = TRUE, layer = "fr")

# Pollution lumineuse
r <- raster("data-raw/World_Atlas_2015.tif")


# Boucle sur les 5 centres CNRS
for (i in 1: nrow(cent)){
  
  # Extraction du centre et géoréferencement
  cnrs <- cent[i,]
  cnrs <- st_as_sf(cnrs, coords = c("longitude", "latitude"), crs = 4326)
  cnrs <- st_transform(cnrs, 2154)

  # Buffer
  cnrsBuf <- st_buffer(cnrs, dist = buf)
  
  #  Population
  comBuf <- st_intersection(com, cnrsBuf) # Intersecter avec le buffer 
  comBuf <- comBuf[,c("GISCO_ID", "LAU_NAME", "POP_2019")] # Ne garder que la population
  countryBuf <- st_intersection(country, cnrsBuf) # Intersecter avec les contours pays
  
  # Raster pollution lumineuse
  rBuf <- st_transform(cnrsBuf, projection(r)) # projection
  rBuf <- intersect(r, rBuf) # intersection raster bbox r30k
  dots <- st_as_sf(as(rBuf, 'SpatialPointsDataFrame')) # transformation vecteur
  dots <- st_transform(dots, crs = 2154) # reprojection
  grid <- st_make_grid(dots, cellsize = 1000, square = TRUE) # grille 1000m
  grid <- aggregate(dots["World_Atlas_2015"], grid, mean) # agrégation valeurs (moy)
  grid <- st_intersection(grid, cnrsBuf) # intersection grille buffer
  
  
  # OSM (Zones d'intérêt touristique)
  bbox <- opq(st_bbox(st_transform(cnrsBuf, 4326))) # Emprise géo de recherche 
  
  # Création d'un dataframe cible
  osm <- data.frame(matrix(ncol = 4, nrow = 0))
  colnames(osm) <- c("osm_id", "key", "value")

  # Exécuter une requête osm et l'associer au dataframe si elle contient 1 objet au min
  osm_q <- function(x, key, value){
    q <- add_osm_feature(opq = bbox, key = key, value = value)
    q <- osmdata_sf(q)
    q <- q$osm_points
    
    if(nrow(q) > 0){
      q$key <- key
      q$value <- value
      x <- rbind(x, q[,c("osm_id", "key", "value")])
    }
    if(nrow(q) == 0){
      x <- x
    }
    return(x)
  }
  
  # 5 requêtes
  osm <- osm_q(x = osm, key = "historic", value = "fort")
  osm <- osm_q(x = osm, key = "historic", value = "monument")
  osm <- osm_q(x = osm, key = "tourism", value = "museum")
  osm <- osm_q(x = osm, key = "tourism", value = "viewpoint")
  osm <- osm_q(x = osm, key = "tourism", value = "artwork")
  
  osm <- st_transform(osm, 2154) # reprojection
  osm <- st_intersection(osm, cnrsBuf)

  # Export couches préparées
  write_sf(obj = cnrs, dsn = paste0("data/", cent[i, "file"], ".gpkg"),
           delete_layer = TRUE, append = TRUE, layer = "cnrs")
  write_sf(obj = cnrsBuf, dsn = paste0("data/", cent[i, "file"], ".gpkg"),
           delete_layer = TRUE, append = TRUE, layer = "cnrsBuf")
  write_sf(obj = comBuf, dsn = paste0("data/", cent[i, "file"], ".gpkg"),
           delete_layer = TRUE, append = TRUE, layer = "com")
  write_sf(obj = countryBuf, dsn = paste0("data/", cent[i, "file"], ".gpkg"),
           delete_layer = TRUE, append = TRUE, layer = "country")
  write_sf(obj = grid, dsn = paste0("data/", cent[i, "file"], ".gpkg"),
           delete_layer = TRUE, append = TRUE, layer = "grid")
  write_sf(obj = osm, dsn = paste0("data/", cent[i, "file"], ".gpkg"),
           delete_layer = TRUE, append = TRUE, layer = "osm")
  
} 

centre <- "oleron"

# Couches disponibles
st_layers(paste0("data/", centre, ".gpkg"))

# Import des couches
cnrs <- st_read(dsn = paste0("data/", centre, ".gpkg"), layer = "cnrs", 
                quiet = TRUE)
buf <- st_read(dsn = paste0("data/", centre, ".gpkg"), layer = "cnrsBuf", 
               quiet = TRUE)
com <- st_read(dsn = paste0("data/", centre, ".gpkg"), layer = "com", 
               quiet = TRUE)
country <- st_read(dsn = paste0("data/", centre, ".gpkg"), layer = "country",
                   quiet = TRUE)
grid <- st_read(dsn = paste0("data/", centre, ".gpkg"), layer = "grid", 
                quiet = TRUE)
osm <- st_read(dsn = paste0("data/", centre, ".gpkg"), layer = "osm",
               quiet = TRUE)

# Ordonner les communes selon leur population (décroissant)
com <- com[order(com$POP_2019, decreasing = TRUE),]

# Carte principale
mf_theme("nevermind") # Un thème 
mf_map(buf, col = "lightblue", border = NA) # Périmètre d'analyse
mf_map(com, col = "white", add = TRUE) # Afficher les communes
mf_map(com, type = "prop", var = "POP_2019", leg_pos = "topleft", 
       col = "grey34", border = "white",
       leg_title = "Population 2019") # Population des communes
mf_map(cnrs, pch = 17, col = "darkred", cex = 2, 
       add = TRUE) # Localisation centre de vacances
mf_label(x = cnrs, var = "centre", halo = TRUE, cex = 0.9, pos = 1, 
         overlap = FALSE, bg = "white", lines = FALSE,
         r = 0.15) # Afficher son nom
mf_label(x = com[1:5,], var = "LAU_NAME", halo = TRUE, cex = 0.8, overlap = FALSE, 
         lines = FALSE, bg = "white") # Label des 5 communes les plus peuplées 
mf_map(buf, col = NA, lwd = 2, add = TRUE) # Périmètre d'analyse

# Carton de localisation
mf_inset_on(fr, pos = "topright", cex = 0.2) 
mf_map(fr, border = NA)
mf_map(fr, col = "white", border = NA, add = TRUE)
mf_map(cnrs, pch = 17, col = "darkred", cex = 2, add = TRUE)
box(col = "white")
mf_inset_off()

# Habillage
mf_layout(title = paste0(round(sum(com$POP_2019),-2),
                         " habitants dans les bourgs des environs"),
          credits = paste0("Ecole Thématique SIG-R, RIATE, 2021\n",
                           "Source : Eurostat - GISCO, 2021"),
          arrow = FALSE)

# Discrétisation et palette de couleur
brks <- c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 1, 1.5, 2, 
          max(grid$World_Atlas_2015, na.rm = TRUE))
pal <- hcl.colors(n = 9, palette = "Lajolla", rev = TRUE)

# Cartographie
mf_map(buf, col = "lightblue", border = NA)

mf_map(grid, var = "World_Atlas_2015", type = "choro", breaks = brks, 
       border = NA, pal = pal, leg_pos = "topleft",
       leg_title = "Pollution lumineuse en 2015\n(Radiance en nano W/cm²)",
       add = TRUE, leg_no_data = "Valeur maximale")

mf_map(grid[which.max(grid$World_Atlas_2015),], col = NA, add = TRUE) # Max
mf_map(country, border = "white", col = NA, add = TRUE) # Front de mer
mf_map(cnrs, add = TRUE, pch = 17, cex = 2, col = "#ffffff90")
mf_map(buf, col = NA, lwd = 2, add = TRUE) # Périmètre d'analyse

# Histogramme (inset)
mf_inset_on(fig = c(0.75, .99, 0.75, 1))
par(mar = c(0,0,1.7,0))
hist(grid$World_Atlas_2015, breaks = brks, col = pal, border = "white", axes = F, 
     labels = "",  xlab="", ylab = "", main ="")
axis(side = 1, at = round(brks, 2), las = 2, tick = FALSE, line = -.9, 
     cex.axis = .7, col.axis = "white")
title("Histogramme de répartition", cex.main = .8, col.main = "white", 
      font.main = 1, adj = 0)
mf_inset_off()

# Mise en page
mf_layout(title = paste0("Voir les étoiles à ", cnrs$lieu, " (Moyenne = ",
                         round(mean(grid$World_Atlas_2015, na.rm = TRUE), 2), ")"),
          credits = paste0("Ecole Thématique SIG-R, RIATE, 2021\n",
                           "Source : Falchi et al., Supplement to the New World",
                           "Atlas of Artificial Night Sky Brightness. GFZ Data Services"), 
          arrow = FALSE)